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ABSTRACT 

The Transit Timing Variation (TTV) method relies on monitoring changes 
in timing of transits of known exoplanets. Non-transiting planets in the system 
can be inferred from TTVs by their gravitational interaction with the transit- 
ing planet. The TTV method is sensitive to low- mass planets that cannot be 
detected by other means. Here we describe a fast algorithm that can be used 
to determine the mass and orbit of the non-transiting planets from the TTV 
data. We apply our code, ttvim.f, to a wide variety of planetary systems to 
test the uniqueness of the TTV inversion problem and its dependence on the 
precision of TTV observations. We find that planetary parameters, including the 
mass and mutual orbital inclination of planets, can be determined from the TTV 
datasets that should become available in near future. Unlike the radial velocity 
technique, the TTV method can therefore be used to characterize the inclination 
distribution of multi-planet systems. 

Subject headings: Stars: planetary systems — celestial mechanics 



1. Introduction 



In Nesvorny & Morbidelli (2008) and Nesvorny (2009) (hereafter NM08 and N09) we 
developed and tested a fast inversion method that can be used to characterize planetary 
systems from the observed Transit Timing Variations (TTVs; Agol et al. 2005; Holman & 
Murray 2005). See NM08 and N09 for a technical description of the method. Here we use 
this new method to solve the TTV inversion problem for an arbitrary planetary system. The 
results provide a baseline for studies of real exoplanetary systems for which TTVs will be 
detected. Examples of past work that would greatly benefit from the application of the fast 
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inversion algorithm discussed here include Steffen & Agol (2005), Agol & Steffen (2007), 
Miller-Ricci et al. (2008) and Gibson et al. (2009). 

In §2, we briefly describe the TTV inversion method. In §3, we apply it to a case 
with coplanar planetary orbits. Inclined planetary orbits are discussed in §4. We show, 
for example, that the mutual inclination of planetary orbits can be determined from TTVs. 
This important parameter, which may be used to test planet-migration theories (e.g., Rasio 
& Ford 1996, Goldreich & Sari 2003), is not typically available from other existing planet- 
detection methods. 



Our TTV inversion method, hereafter TTVIM, has two parts. The first part is a fast 
algorithm for the computation of transit times, (Stj) tT iai, 1 < J ' < N, for specified planetary 
parameters. This algorithm is based on perturbation theory (NM08, N09). It calculates 
the short-period TTVs as these have been shown to be the most diagnostic (NM08). The 
long-term effects such as the apsidal precession produced by the perturbing planet are more 
difficult to detect if transit observations span only a few years (Miralda-Escude 2002; Heyl 
& Gladman 2007). 

The second part of TTVIM is an adaptation of the Downhill Simplex Method (DSM; 
Press et al. 1992). The DSM is used to search for the minima of 



where (<5t,)t r iai are the transit times produced by the first part of the algorithm for a trial 
planetary system, (<5t,) bs are the observed mid-transit times, and <jj are the measurement 
errors. It is assumed here (as indicated by 5's) that the period of the transiting planet, P, 
has been removed from transit observations. Thus, (Stj) ohs = (tj) ohs — CjP, where (t,) bs are 
the actual transit times and integer Cj denotes the transit cycle. 

The best-fit planetary parameters correspond to the global minimum of x 2 , denoted by 
Xmin in the following. A large number of initial trials must be used to assure that the DSM 
method finds Xmin- The confidence levels for the normally distributed data can be defined 
as Ax 2 = X 2 — Xmin < (Ax 2 )cut, where the (A% 2 )cut values are properly chosen for N and 
the required confidence level (NM08). 

Here we assume that the mass and orbit of the transiting planet are known from transit 
and radial velocity (RV) measurements as this should be the most common case in practice. 



2. 



Inversion Method 




(1) 
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If so, x 2 is a function of seven unknown parameters of the perturbing planet, 
X 2 — X 2 ( m 2, 0,2, e 2 , % 2 , ^2, A 2 ), where m 2 is the mass, a 2 semimajor axis, e 2 eccentricity, 
i 2 inclination, Q 2 nodal longitude, w 2 periapse longitude, and A 2 the mean orbital phase at 
t — (arbitrarily defined here to correspond to cycle C*o)Jj The parameters of the transiting 
planet will be denoted by index 1. DSM must therefore search in 7D space for the global 
minimum of \ 2 ■ This is not a trivial task because \ 2 often has many deep and narrow local 
minima (Steffen & Agol 2007). Fortunately, several simplifications can be made. 

First, as the amplitude of the short-period TTVs scales linearly with m 2 , we can calcu- 
late the TTV profile for the selected a 2 ,e 2 ,i 2 ,Q 2 ,zu 2 , X 2 values and obtain m 2 by the linear 
least-square fit. Second, the determination of x 2 f° r a new set of the Q 2 , w 2 and A 2 values 
is computationally cheap in the perturbation algorithm, if \ 2 was determined previously for 
the required a 2 , e 2 and i 2 valuesEI The code can thus efficiently search for the minimum of 
X 2 (o 2 , e 2 , i 2 ) for any value of Q 2 , zu 2 and A 2 . In practice, we use 5 to 20 values between and 
360° to resolve each of these parameters. This effectively reduces the number of dimensions 
to 3. Once the interval of estimated a 2 , e 2 and i 2 values is narrowed down, the solution can 
be refined by using the full DSM search in 7D. 

The tricky part of TTVIM is the choice of the initial guess in the (a 2 ,e 2 ,i 2 ) space. 
By trials and errors we found that fine (and non-linear) sampling of a 2 is generally needed 
for a successful convergence of the algorithm. The best results were obtained with uniform 
sampling in I /a 2 , where a — a\/a 2 < 1. Parameters e 2 and i 2 require less care since DSM 
usually finds the right minimum even in the high-e 2 and/or high-z 2 case if at least one corner 
of the initial simplex is stretched to e 2 > 0.2 and i 2 > 30°. 

With the nominal setup, our TTVIM code (ttvim.f) requires about 2 minutes of CPU 
tim^l for a coplanar fit with i 2 = and about 50 minutes for the full 7D fit. In the absence 
of measurement errors, the success rate in finding xLm is better than 95%. Thus, ttvim.f is 
a robust code that can reliably solve the TTV inverse problem at a low computational cost. 



1 These are the actual parameters used in DSM. The boundary at e2 = does not need a special treatment 
because (e2,W2) is formally equivalent to (— e2,vj2 +n/2). Similar rules apply to (22,^2) and (m2,A2). 

2 This is because all Fourier terms can be pre-computed for et 2 , & n d ii and need only to be assembled 
with the specific O2, ^2 and A2 values. The assembling procedure itself is computationally inexpensive. 

3 On a single 2.7-GHz Opteron processor. 



-4- 



3. Results for Coplanar Orbits 

We used a random number generator to define different sets of parameters m 2 , a 2 , e 2 , 
i 2 , ©2 and A2. Typically, f 000-2000 different planet parameter sets were used in tests. 
In each case, the orbital evolution of the two planets was followed for a fixed timespan, 
< t < T int , with the Bulirsch-Stoer integrator (Press et al. 1992). During this timespan we 
interpolated for and recorded all transit times of the inner planet. This data mimics the real 
observations, (<%) bs- They were used in a blind test where we applied the TTVIM code to 
each of these cases in an attempt to recover the original mass and the orbital elements of 
the non-transiting planet. 

We start by discussing the case with star's mass m = Ms U n, where Ms un = 2 x 10 33 g 
is the mass of the Sun, mi = 10~ 3 mo, a\ = 0.1 AU, e\ — i\ — and N = 100 consecutive 
transits. Since NM08 and N09 showed that the behavior of the inversion method is insensitive 
to mi, we will not test different mi values in this work. To distinguish between the issues 
related to the intrinsic limitations of ttvim.f and those arising from the finite precision of 
the real measurements, we first discuss an idealized case with zero measurement noise. 

For coplanar orbits and Oj = 0, the TTVIM code finds the correct planetary parameters 
with a high rate of success (Fig. [Q. The typical precision in the successful cases is |m 2 — 
m* 2 \/m 2 < 0.2, \a 2 -a* 2 \/a 2 < 0.02, |e 2 — e* 2 \ < 0.02, \w 2 -wl\ < 10° and | A 2 — A| | < 10°, where 
the asterisk denotes the values determined by the TTVIM codeB This is very satisfactory. 
In the absence of measurement errors, the result of the TTVIM code illustrated in Fig. [1] 
with m 2 = 10~ 4 mo is insensitive to the actual value of m 2 . 

The main failure mode of the TTVIM code occurs near mean motion resonances between 
planets, because resonant perturbations are not (yet) taken into account in TTVIM. While 
the resonant signal can improve our chances of the TTV detection for (near-) resonant planets, 
it seems less useful in helping us estimate the mass and orbit of the planetary companion. 
Specifically, the amplitude and period of the resonant signal can be fit by a number of 
different planetary setups corresponding to different resonances. Thus, without an apriori 
knowledge of the mean motion resonance that is responsible for the observed behavior, the 
inversion problem from resonant frequencies alone is strongly degenerate. 

Fortunately, the short-period TTVs underlying the resonant signal can still be used to 
determine the planetary parameters without much ambiguity. As shown in NM08, probably 
the best strategy is to isolate short-period frequencies in the signal by Fourier filter and 
apply the inversion method to the filtered signal. The application of this procedure is 



4 Except for very small values of e-i for which the errors in W2 can be large. 
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straightforward in individual cases (see NM08), where the resonant period, and thus the 
appropriate frequency cutoff, can be estimated from (t,-) bs- We verified that this procedure 
works quite well in >75% of cases shown in Fig. [T]in which the resonant variations are an 
issue. 

The remaining <25% unsuccessful cases (representing <5% overall) correspond to the 
very large values of e 2 for which the Laplacian expansion of the perturbing function in 
TTVIM is not convergent (NM08), and/or planetary configurations that are not Hill stable. 
Direct iV-body integrations can be used to address the TTV inversion problem in the very- 
high-eccentricity domain, but the CPU cost of these tests is likely to be substantial and lies 
beyond the scope of this letter. 

The measurement errors have a profound effect on the uniqueness of the inverse problem 
(Fig. [2]). For m 2 = 10~ 4 m , N = 100 and Oj = a = 3 s, corresponding to the Kepler- 
like precision of timing measurements for a Sun-like star with a 2 Neptune-mass planet, 
unique determination of planetary parameters can be achieved for most stable systems with 
<?2 = fl 2(l — e 2 ) < 3.3ai, while for aj = a = 10 s (Corot-like precision), it is required that 
g 2 < 2.6a!. These limits approximately correspond to the planetary parameters for which 
the amplitude of the short-period TTVs is comparable to a. 

Figure [3] shows the result of the TTVIM code for an Earth- mass planet. The region 
of parameter space in which unique determination can be achieved from TTVs is relatively 
small even with a = 1 s. Thus, an Earth-mass planet detection and characterization of its 
orbit will require a rather fortuitous setup of the planetary system, in which (g 2 — < 2 

(for external perturber). 

4. Results for Inclined Planetary Orbits 

We applied the TTVIM code to mock planetary systems with < z 2 < 50°. As in §3, 
we assumed that m = M Sun , mi = 10~ 3 m , a± = 0.1 AU, e x = and used iV = 100 
consecutive transits. Figure 0] shows the result for a,j = a = 0. The (a, e) plot does not 
differ much from the coplanar case although it may be noted that the quality of fits slightly 
degraded for the large e 2 values. This is probably related to the convergence problems of 
the perturbation algorithm in TTVIM. A precise iV-body integrator should perform better 
for high e 2 , although it has yet to be shown that an iV-body integrator can be applied to 
the inclined inverse problem in practice due to the large CPU cost. 

Probably the most exciting result obtained in this work is that it was possible to de- 
termine the mutual inclination of planets for most planetary systems (Fig. Hb). Unlike the 
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radial velocity technique, the TTV method can therefore be used to characterize the incli- 
nation distribution of multi-planet systems. Figure [5] shows the detailed statistic of TTVIM 
errors in i 2 . In most cases, |« 2 — £j>| < 2°. The tail of larger \i 2 — i 2 \ values corresponds to 
the high-eccentricity cases. If the statistic is limited to q 2 = a 2 (l — e 2 ) > 0.25 AU (Fig. Eh; 
dashed line), the fraction of successful cases with \i 2 — i 2 \ < 2° increases to >90%. In the 
succesful cases, orbital angles Q 2 , ^2 and X 2 are generally correctly determined to within a 
better than 5° precision. 

We also studied how the uniqueness of the inclined inverse problem is affected observa- 
tional errors. The trends seen in these tests are very similar to those described in §3. Namely, 
the instrumental noise sets an upper limit on q 2 beyond which the determination of plan- 
etary parameters from TTVs is ambiguous. Again, we see that these limits approximately 
correspond to the planetary parameters for which the amplitude of the short-period TTVs 
is comparable to a. The results in N08 and Payne et al. (2009) can therefore be used to 
estimate whether (or not) a unique characterization of the specific inclined planetary system 
may be achieved from TTV observations with given a. 

We find that TTVs obtained in the coplanar case represent a good approximation of the 
TTVs for planetary orbits with i 2 < 20°. The planar version of the TTVIM algorithm can 
therefore be used in these cases to estimate the a 2 and e 2 values of the perturbing planet. 
This helps to narrow the range of initial guesses for the 7D fit and represents a factor of ~20 
speed up of the inversion. The full 7D algorithm needs to be used for i 2 > 20°. 

5. Conclusions 

The method developed here can be used to analyze TTVs found for any of the potentially 
hundreds of planets expected to be discovered by Kepler (Beatty & Gaudi 2008). Kepler 
should be able to detect transit timing variations of only a few seconds (Holman & Murray 
2005), which should easily exist in many systems, extrapolating from the radial velocity 
planets (Agol et al. 2005, Fabrycky 2008). 

Perhaps the most interesting result that comes out of this work is that the shape of the 
TTV signal is generally sensitive to the orbital inclination of the non-transiting planetary 
companion. Thus, the TTV method can provide means of determining mutual inclinations 
in systems in which at least one planet is transiting. This parameter cannot be determined 
by other planet-detection methods. 

TTVIM algorithm can be easily extended to incorporate uncertainties in the transiting 
planet's parameters. This can be done by sampling dimensions that correspond to the 
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additional parameters. For example, in N09 we extended the NM08 method to the case with 
e\ 7^ 0. This may be especially relevant to the transiting planets that will be discovered by 
Kepler because these planets are expected to have wider orbits, which are less susceptible 
to the circularizing effects of tides. The low CPU cost of the TTVIM algorithm is the key 
element which will make such studies possible. 

This work was supported by the NSF AAG program. 

REFERENCES 

Agol, E., Steffen, J. 2007, MNRAS, 374, 941 

Agol, E., Steffen, J., Sari, R., Clarkson, W. 2005, MNRAS, 359, 567 

Beatty, T. C, Gaudi, B. S. 2008, ApJ, 686, 1302 

Fabrycky, D. C. 2009, IAU Symposium, 253, 173 

Gibson, N. P., et al. 2009, MNRAS, 1684 

Goldreich, P., Sari, R. 2003, ApJ, 585, 1024 

Heyl, J. S., Gladman, B. J. 2007, MNRAS, 377, 1511 

Holman, M. J., Murray, N. W. 2005, Science, 307, 1288 

Miller-Ricci, E., et al. 2008, ApJ, 682, 586 

Miralda-Escude, J. 2002, ApJ, 564, 1019 

Nesvorny, D., Morbidelli, A. 2008 (NM08), ApJ, 688, 636 

Nesvorny, D. 2009 (N09), ApJ, 701, 1116 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. 1992, Numerical recipes 
in FORTRAN. The art of scientific computing. (Cambridge: Cambridge University 
Press) 

Rasio, F. A., Ford, E. B. 1996, Science, 274, 954 
Steffen, J., Agol, E. 2005, MNRAS, 364, L96 



- 8- 

Steffen, J. H., & Agol, E. 2007, Transiting Extrasolar Planets Workshop, 366, 158 
fla7Xiv:astro-ph/0612442p 



This preprint was prepared with the AAS IATgX macros v5.2. 



- 9- 




Semimajor axis (AU) 



Fig. 1. — TTVIM code results for planetary systems with m = Ms un , m\ = 10 _3 mo, 
a\ = 0.1 AU, e\ = and z 2 = 0. Planetary parameters for which the TTVIM code converged 
to the correct solution were denoted by blue x's. Incorrect solutions were denoted by red 
dots. We defined the correct solution as having |m 2 — m 2 |/m 2 < 0.5, |a 2 — a 2 |/a 2 < 0.05 and 
\ e 2~ 62 1 < 0.05, where m 2 , a 2 and e 2 are the original planetary parameters for which the TTV 
signal was computed by iV-body integration, and m 2 , a 2 and e 2 are the values determined by 
the TTVIM code. In the majority of cases corresponding to correct solutions, the TTVIM 
code determined the original orbital parameters with a better than 2% precision and mass 
with a better than 20% precision. The two bold solid lines show the planet- crossing (upper) 
and Hill- stability limits (lower). We also show the location of the principal mean motion 
resonances between the planets (e.g., 2:1 at a 2 = 0.16 AU). There are two lines per resonance 
corresponding to the left and right separatrices of resonant motion. The V-shaped profiles 
are characteristic for mean motion resonances that become wider with eccentricity. 
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Fig. 2. — TTVIM code results for a planet with m-i = 10 4 mo and different levels of the 
measurement error, a. See caption of Fig. [1] for a description of different lines and symbols. 
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Fig. 3. — TTVIM code results for an Earth- 
different levels of the measurement error, a. 
different lines and symbols. With o > 3 s, 
Earth-mass planets with very specific orbits. 



mass planet [m-i = 3 x 10 -6 Ms un ) and two 
See caption of Fig. CD for a description of 
the TTVIM code can only characterize the 
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Semimajor axis (AU) 

Fig. 4. — The same as Fig. [1] but with i 2 ^ 0. Most correct solutions (blue x's) have 
\a 2 - a* 2 \/a 2 < 0.015, |e 2 - e* 2 \ < 0.02, \i 2 - i 2 \ < 2° and \m 2 - m* 2 \/m 2 < 0.2. In (b), the 
dashed lines denote the libration centers of the principal mean motion resonances. 
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Inclination Error (deg) Relative Mass Error 



Fig. 5. — Distribution of TTVIM errors in % 2 (left) and m 2 (right) for the case shown in 
Fig. HI We show the total distribution (solid line) and the one for q 2 > 0.25 AU (dashed). 
In the later case, the erroneous determinations with \i 2 — i 2 \ > 5° are reduced because the 
algorithm does not need to deal with the difficult case when q 2 ~ ai. Most cases correspond 
to |m 2 — ra1^\jm2 < 0.2 (i.e., <20% precision of mass determination) and |z 2 — i 2 \ < 2°. 



